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Abstract. We describe and demonstrate a method to reconstruct an amplitude 
*"^ , equation from the nonlinear relaxation dynamics in the succession of the Rosensweig 

instability. A flat layer of a ferrofluid is cooled such that the liquid has a relatively 
high viscosity. Consequently, the dynamics of the formation of the Rosensweig pattern 
becomes very slow. By sudden switching of the magnetic induction, the system is 
pushed to an arbitrary point in the phase space spanned by the pattern amplitude 
and the magnetic induction. Afterwards, it is allowed to relax to its equilibrium point. 
■^j- ' From the dynamics of this relaxation, we reconstruct the underlying fully nonlinear 

equation of motion of the pattern amplitude. The measured nonlinear dynamics serves 
to select the best weakly nonlinear expansion which describes this hysteretic transition. 
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1. Introduction 

When Aristotle (350 B.C.) described the motion of falling bodies, he incorrectly claimed 
that the velocity is constant and proportional to the mass of the body, because he did 
not know the concept of inertia. He has been later proven wrong by Galilei (1638) 
and Newton (1687), who linked acceleration with forces. Aristotles' idea of motion, 
however, is a valid approximation in many cases, when a viscous fluid is involved and 
either the typical dimension d and the velocity v are small, or the kinematic viscosity 
v is high. In this limiting case of small Reynolds numbers Re = vdjv <C 1, the viscous 
forces outweigh inertia, and body and fluid motion is determined entirely by drag forces. 
Therefore, from a measurement of the velocity field one immediately obtains the acting 
forces. 

A similar dynamics is expected in the neighbourhood of a hydrodynamic instability, 
namely an amplitude equation of the form dA/dt = f(A). The coefficients of these 
functions / can be determined experimentally. As an early example, Wesfreid et al. 
(1978) accomplished this in a Rayleigh-Benard convection experiment by observing the 
relaxation of the pattern amplitude following a jump in the control parameter. In this 
way, the coefficients of the amplitude equation of a pitchfork forward bifurcation to 
stripe-like patterns had been obtained. Here we present such a measurement technique 
for a transcritical bifurcation to hexagons, which is characteristic for systems with broken 
up-down symmetry, i.e. a different universality class. As a specific example, we study 
the hexagonal arrangement of spikes forming in a layer of magnetic fluid. 

Magnetic fluids, also known as ferrofluids, are colloidal suspensions of magnetic 
nanoparticles (Rosensweig 1985). When a vertically oriented magnetic field is applied 
to a horizontal layer of ferrofluid, then a regular pattern of liquid spikes can appear above 
a certain threshold value of the magnetic induction B c . This normal field instability was 
first reported by Cowley & Rosensweig (1967) and is thus also known as the Rosensweig 
instability. For an infinitely deep container, these authors provide also a linear stability 
analysis to find the critical threshold B c of the magnetic induction and the critical wave 
number k c . Their analysis works by computing the dynamics of a regular pattern in the 
vicinity of the threshold induction and for very small amplitudes. This approach has 
later been extended to calculate the viscosity dependent growth rate by Salin (1993), 
and a finite depth of the container was taken into account by Weilepp & Brand (1996). 
For the case of a viscous magnetic fluid and an arbitrary layer thickness, Lange et al. 
(2000) derived the growth rate and the wave number of maximal growth. The latter is 
experimentally determined and compared with the theoretical predictions in the same 
article, and later on compared with the wave number of the nonlinear state (Lange 
et al. 2001). Lange (2001) calculates the growth rate for ferrofluids of different viscosity 
and eventually extends his analysis to the case of a nonlinear magnetization curve 
(Knieling et al. 2007). 

The nonlinear aspects of the dynamics have posed additional difficulties to the 
theoretical analysis (Lange 2002). Only recently, Bohlius et al. (2007, 2008) successfully 
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tackled the issue of systematically deriving nonlinear amplitude equations. For the 
amplitude A of hexagonal patterns these authors obtained 

T ^ = eA + >y 1 A*-gA*. (1) 

Here, To denotes a time scale, e = (B 2 — B 2 )/B 2 measures the distance from the 

bifurcation point, and 71 and g are expansion coefficients. The static solution of this 

amplitude equation coincides in structure with the result presented by Gailitis (1977) 

and by Kuznetsov & Spektor (1976), which was derived for magnetic susceptibilities 

X < 1 using an energy minimization approach. Using a similar approach, Friedrichs 

& Engel (2001) present a solution which is applicable up to x < 0.4. This solution is 

slightly different and corresponds to the amplitude equation with 

r)A 
To — = eA + ll (l + e)A 2 ~gA 3 . (2) 

The difference to the cubic equation (JTJ) is the additional factor (1 -+- e) in the quadratic 
term. The physical reasoning behind this higher order correction is a scaling of the 
quadratic coefficient with B 2 , which is proportional to 1 + e = B 2 /B 2 . In the limit 
of s — > 0, both equations coincide. Now the question arises, which of these solutions 
compares best to the experiment. 

So far, the static solution of equation fl2]) was successfully fitted to the measured 
pattern amplitude for a ferrofluid with high susceptibility (Richter & Barashenkov 
2005) and moderate susceptibility (Gollwitzer et al. 2007). However, the amplitude 
equation (CQ) provides much more information, specifically the nonlinear dynamics of 
the pattern formation. This has not yet been exploited. Knieling et al. (2007) observed 
the nonlinear dynamics in a previous experiment and in related numerical simulations, 
but a quantitative comparison to an analytical model is performed only in the linear 
regime of the pattern growth. 

The key idea in this paper is now to use a fluid with a very high viscosity to 
study the pattern formation of the Rosensweig instability. Consequently, the Reynolds 
number is very small Re ~ 10~ 3 , and the relaxation is monotonic and slow. The heavily 
overdamped system makes it easy to study the nonlinear dynamics quantitatively: due 
to the strong damping, the pattern amplitude follows a first order equation of motion 

f = /(M) (3) 

even in the nonlinear regime. Here / is an analytic function, which does not directly 
depend on the time t. From the experimental data, / can be reconstructed and directly 
compared to the theoretical amplitude equations. 

In the next section, we describe the experimental setup and the magnetic fluid we 
use. Subsequently, the measurement protocol and the data analysis are described. Then 
the experimental data are compared to theoretical models, and finally the results are 
discussed. 
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Figure 1. Setup of the apparatus for dynamic measurements of the Rosensweig 
instability. Scheme of the assembled setup (a), and exploded view of the container 
with the coils (b). 



2. Experimental setup 



The experimental setup for the measurements of the surface deflections consists of an 
X-ray apparatus described in detail by Richter & Biasing (2001). An X-ray point source 
emits radiation vertically from above through the container with the ferrofluid, which is 
placed midway between a Helmholtz pair of coils. Underneath the container, an X-ray 
camera records the radiation passing through the layer of ferrofluid. The intensity at 
each pixel of the detector is directly related to the height of the fluid above that pixel. 
Therefore, the full surface topography can be reconstructed after calibration (Gollwitzer 
et al. 2007). 

The container, which holds the ferrofluid sample, is depicted in figure ^b). It 
is a regular octagon machined from aluminium with a side length of 77 mm and two 
concentric inner bores with a diameter of 140 mm. These circular holes are carved from 
above and below, leaving only a thin base in the middle of the vessel with a thickness of 
2 mm. On top of the octagon, a circular lid is placed, which closes the hole from above 
(see figure []]&). Each side of the octagon is equipped with a thermoelectric element 
QC-127-1.4-8.5MS from Quick-Ohm. These are powered by a 1.2 kW Kepco KLP-20- 
120 power supply. The hot side of the thermoelectric elements is connected to water 
cooled heat exchangers. The temperature is measured at the bottom of the aluminium 
container with a PtlOO resistor. A closed loop control, realized using a computer and 
programmable interface devices, holds the temperature 9 of the vessel constant with an 
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Table 1. Material properties of the ferrofluid mark APGE32 (Lot G090707A) from 
Ferrotec Co. 



Quantity 



Value 



Error 



Unit 



Density 
Surface tension 
Viscosity at 10 °C 
Saturation magnetization 
Initial susceptibility at 10 °C 



P 

(7 

n 

M s 

X.Q 



1168 
30.9 

4.48 
26.6 



±1 
±5 
±0.1 
±0.8 



3.74 ±0.005 



kgm 3 
mNm -1 
Pas 
kAm" 1 



Fit of M(H) with the model by Ivanov & Kuznetsova (2001): 
Exponent of the r-distribution 
Typical diameter 
Volume fraction 



a 
do 



3.8 

1.7 

5.96 



±1 

±0.2 

±0.2 



critical induction from M{H) 



B r 



10.5 ±0.1 



mT 



accuracy of 10 mK. 

The container is surrounded by a Helmholtz pair of coils, thermally isolated from 
the vessel with a ring made from a flame resistant composite material (FR-2). The 
size of the coils is adapted to the size of the vessel in order to introduce a magnetic 
ramp. With these coils, the magnetic field strength falls off towards the border of the 
vessel, where it reaches 80 % of its value in the centre. For details of the magnetic ramp 
see Knieling et al. (2010). Filling the container with ferrofluid enhances the magnetic 
induction in comparison with the empty coils for the same current /. Therefore B(I) 
is measured immediately beneath the bottom of the container, at the central position, 
and serves as the control parameter in the following. 

We use the commercial ferrofluid APG E32 from Ferrotec Co, the basic material 
parameters of which are listed in table [TJ The density was measured using a DMA 
4100 density meter from Anton-Paar. More problematic is the measurement of the 
surface tension. We measure the surface tension using a commercial ring tensiometer 
(LAUDA TE 1) and a pendant drop method (Dataphysics OCA 20). Both methods 
result in a surface tension of a = 31 ± 0.5mN/m, but when the liquid is allowed to 
rest for one day, a drops down to 25 ± 0.5mN/m. This effect, which is not observed in 
similar, but less viscous magnetic liquids like APG 512 a, which was used in previous 
work (Gollwitzer et al. 2009), gives a hint that our liquid is chemically less stable in 
the sense that the surfactants change the surface tension on a longer time scale, when 
the surface is changed. Since the pattern formation experiments do change the surface 
during the measurements, the uncertainty of the surface tension is about 5 mN/m, as 
given in table [TJ 

The viscosity rj deserves a special attention for the experiments in this paper, as it 
influences the time scale of the pattern formation. It has been measured in a temperature 
range of — 5°C < 9 < 20 °C using a commercial rheometer (MCR301, Anton Paar) in 
cone and plate geometry (figure [2]). At room temperature, the ferrofluid we use has a 
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Figure 2. The dynamic viscosity 
r\ of the ferrofluid APG E32 versus 
the temperature 0. The symbols 
represent measurements, and the 
solid line is an approximation by 
equation (J4j) - 
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Figure 3. Magnetization curve 
of the ferrofluid APGE32. The 
symbols show the measured data 
at 6 = 20 °C. The black dashed 
line marks a fit with the model 
by Ivanov & Kuznetsova (2001). 
The blue solid line shows an 
extrapolation to 9 = 10 °C 
according to this model. 



viscosity 7/ = 2 Pas, i.e. 2000 times the viscosity of water. This value can be increased 
by a factor of 9 when the liquid is cooled to —5 °C. The data from figure [2] can very well 
be fitted with the Vogel-Fulcher law (Rault 2000) 

r] = r] exp (— —J , (4) 

with 7/0 = 0.48 mPas,j/; = 1074 K, and 9 = — 107.5 °C, as marked by the solid line 
in figure [2j This means, that the viscosity can easily be adjusted in a wide range by 
controlling the temperature. For the present measurements, we chose a temperature of 
9 = 10 °C, where the viscosity amounts to 7/ = 4.48 Pas according to equation @. 

The magnetization curve has been determined using a fluxmetric magnetometer 
consisting of a Helmholtz pair of sensing coils with 6800 windings and a commercial 
integrator (Lakeshore Fluxmeter 480). The sample is held in a spherical cavity with a 
diameter of 12.4 mm in order to provide a homogeneous magnetic field inside the sample 
with a demagnetization factor of |. The magnetization curve has been measured at a 
temperature of 9 = 20 °C (see figure [3]). For a comparison with the pattern formation 
experiments, this curve must be extrapolated to 9 = 10 °C, which is done using the 
model by Ivanov & Kuznetsova (2001). Following Rosensweig (1985), §7.1, from this 
curve and the material parameters a critical induction of 10.5 mT is estimated (c.f. table 

ED. 
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Figure 4. The multistep measurement protocol for the relaxation measurements. 
Dotted arrows indicate jumps of the magnetic inductions. Blue arrows denote the path 
of the system during the relaxation phases. A movie which shows the sequence l-2a 



is available from http://ep5.uni-bayreuth.de/ff/movies/Frozensweig2a.html 



3. Measurement protocol 

Figure H] displays the measurement protocol on the basis of the bifurcation diagram. 
The static pattern amplitude of the Rosensweig instability in our fluid is indicated by 
the red line. When the system is set onto an arbitrary initial point (B ini , Am) in this 
diagram, and the magnetic induction B is kept constant, the amplitude A monotonically 
increases or decreases, until the system reaches the stable equilibrium (solid red line). 
The direction of the change of A depends on the region, where (-Bj m , Am) is situated - 
in the regions I and III in figure HI A increases, and in regions II and IV, the amplitude 
decreases with time. 

In order to push the system to an arbitrary initial location (B ini , Ani), a multistep 
measurement protocol is employed. The first step (path 1 ) is always a relaxation of the 
pattern in region I at the overcritical induction -Bhigh = 11.45 mT for 1 min, to reach the 
final amplitude of Aiigh = 2.98 mm at that point. The corresponding pattern is shown 
in figure [H a). To get to an arbitrary point in the region II in diagram HI a second 
step is needed. When the system is settled at (-B m gh 5 ^mgh) 5 the magnetic induction is 
quickly changed to the desired value B ini , and the resulting dynamics is observed (path 
2a), until the final amplitude is reached. To get into the inner regions III or IV, three 
steps are needed in total - coming from (-Bhigh, Aiigh), the second step (path 2b) is a 
decay of the pattern amplitude at the subcritical induction B\ ow = 10.74 mT, until the 
intended initial amplitude Am is reached. The induction is then quickly raised to the 
desired B m [. The system then follows the path 3a or 3b in region III or IV, respectively. 

When the desired initial amplitude Ani is zero, for example when the growth of 
the pattern in region I shall be observed, we also apply the three step protocol with the 
detour by paths 1 and 2b. At this point, the pattern decays until the amplitude is very 
small (Am = 0.05mm). After that, the induction is raised to B — B ini . We use this 
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Figure 5. The final pattern at B — 11.45 mT. Chart (a) displays a reconstruction of 
the surface in real space. The outer dimension of the container is not to scale. The 
colour code gives the height of the liquid surface above ground in mm. The pattern 
amplitude is determined from the corresponding power spectrum shown in (6) by the 
total power in the encircled mode. 



procedure instead of directly switching the magnetic induction from zero to B = B m [ in 
order to establish the exact same pattern in all regions. Coming from a perfectly flat 
surface, the pattern would have additional degrees of freedom manifesting themselves 
in point defects or different orientations of the wave vectors. When we take the detour 
by the paths 1 and 2b, we seed the wave vectors of the pattern at (-Bhigh, Aiigh), and the 
final pattern developed is likely to be the same. 

We explored all regions in figure H] in several rounds. In the first round, we examined 
region I. We used the three step procedure to observe the growth of the pattern from a 
very small amplitude up to the stable solution for 21 different inductions in the range 
11.19 mT < B < 11.46mT. In the second round we explored region II. Starting from 
Ani = ^4high and 65 different inductions in the range 10.61 mT < B < 11.46 mT, the 
decay of the pattern was observed until it reaches the stable solution or zero. The third 
to sixth rounds covered the regions III and IV. Starting from four different amplitudes, 
we observed the evolution of the pattern in the range of 10.89 mT < B < 11.27mT. 
At each round, we recorded the complete evolution of the surface topography of the 
ferrofluid during the last step with the fastest possible frame rate (7.5 Hz) of the X-ray 
device, in total taking about 170 000 frames. 



4. Data analysis 

4-1. Reconstruction of the pattern amplitude 

After processing the image data, we finally arrive at the surface topography for every 
frame. Figure [5(a) shows a reconstruction of the surface at B = -Bhigh- A description of 
the technical details of this process can be found in the paper by Gollwitzer et al. (2007). 
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Figure 6. The evolution of the pattern amplitude at an induction B = 11.06 mT, 
starting from different amplitudes (a). The labels of the datasets correspond to the 
paths in figure @] The path 3c ends in a localized state with only one spike (X- 
ray image shown in 6), and has been removed from the data. A movie showing 
the temporal evolution of the pattern to the final state in (6) is available from 
http : //ep5 . uni -bayreuth . de/f f /movies/Frozensweig3c . html 



The amplitude of the pattern, A, is then determined in Fourier space (figure [5]&). We 
use a circularly symmetric Hamming window with the weight function 

. ..wg+z^V , 2 L ^.2 

w{x,y) 



x 



y <ri 



(5) 



else 



for apodization (Gollwitzer et al. 2006) with a radius r w = 46 mm. The total power 
in one of the modes (marked with a red circle in figure [5]6) is used to compute the 
amplitude of the pattern 



.4 




(6) 



where Cj are the Fourier coefficients inside the circle. In order to get a meaningful 
estimate, the normalization factor N is chosen such that A is the height difference 
between maxima and minima, when the input is a perfectly sinusoidal hexagonal pattern. 
Figure \6^ a) displays the resulting amplitude versus time for the paths labelled 2a, 
3a and 3b in figure HI As expected from the bifurcation diagram, the amplitudes for 
paths 2a and 3b decrease monotonically, and the stable solution is approached in an 
asymptotic fashion. Similarly, 3a increases monotonically and asymptotically converges 
to the patterned state. There is, however, a small glitch. The bifurcation diagram 0] 
provides only a single solution for the patterned state, but the asymptotic amplitudes of 
the paths 2a and 3a differ by 2.6 %. The reason is that the final pattern is different, in 
spite of the aforementioned efforts - path 2a ends in a pattern with 20 spikes, whereas 
the final pattern of path 3a contains only 11 spikes. These additional spikes are situated 
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Figure 7. Reconstruction of the amplitude equation at an induction B = 11.05 mT. 
The symbols are central differences from equation (|7J| . The red line serves as a guide to 
the eye. The blue arrows indicate the range of the corresponding data sets presented 
inFig.[l](a). 



near the border of the window function and therefore contribute only little to the overall 
amplitude. A further path, 3c, which starts very close from the initial amplitude of the 
path 3b, ends in a finite intermediate amplitude. Path 3c finally ends in a localized 
solution with only one single spike as depicted in figure [^ b) . 

We are currently unable to control these additional degrees of freedom for the 
system. In the following, the datasets with a reduced number of spikes have not been 
treated in a special way. Only when the final pattern does not fill the full width at 
half maximum of the window function, that is for less than 10 spikes, the estimated 
amplitude would be more than 9 % too small for the given window function. These 
datasets have been sorted out. 

Out of 205 datasets, 32 have been rejected, the final pattern of which consists 
of only a small number of spikes. Every pattern from a single spike up to a cluster 
of 9 has been found, whereas the fully developed pattern contains up to 27 spikes. 
These patterns emerged in the bistable regime of the fully developed pattern and have 
been previously observed by Richter & Barashenkov (2005). In their experiment, the 
localized spikes have been initiated by a local disturbance of the magnetic induction. 
In the present experiment, these patterns emerge spontaneously from a starting point 
(.Bini; Am) near the unstable solution branch. In this study, however, we focus only on 
the fully developed patterns. 



4-2. Recovery of the amplitude equation 

In the following, we will describe a method to extract the amplitude equation from the 
experimental data. So far, we have measured the evolution of the amplitude A(t) for 
different values of the induction B and a set of initial amplitudes A[ n \. Suppose, that the 
system can indeed be described by an amplitude equation of the form ([3]). Then A(t) is 
the solution of this equation for different initial conditions (_B ini , y4 ini ), and we want to 
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Figure 8. The measured amplitude dynamics for the whole phase space. The colour 
indicates the velocity dA/dt in mm/s. Red (blue) means rising (falling) amplitude, 
respectively. Green indicates the zero (root of the amplitude equation) . 

get the function / on the right hand side of (J3J) . Because this function should not depend 
on time, the time derivative A = dA/dt of our measured A(t) directly gives the value of 
this function at the corresponding amplitude. We therefore plot A as a function of A in 
figure [7J for one selected induction in the bistable region. The time derivative has been 
estimated from the measured data by central differences (Abramowitz & Stegun 1965, 
§25.1.2) 



dA 
~dt 



n+i 



A\ 



An+l — A n 

At 
A n -\-i + A n 



n+4 



(?) 



where At = 0.134 ms is the time between consecutive frames, and Ai denotes the pattern 
amplitude of the i th frame. The plot in figure [7J comprises three different data sets. 
The sets 2a and 3a stem from a relaxation to a pattern amplitude of about 2 mm; 
set 3b describes the decay to a flat surface layer. Because the time derivative of the 
experimental data suffers heavily from noise, we also plot a smooth approximation 
explained later in § [5] as a guide to the eye (red line in figure [7]). In the bistable 
regime, the amplitude equation has three roots, corresponding to one unstable and two 
stable solutions. In figure [7J these stable and unstable solutions are characterized by 
zero-crossings with a negative or positive slope, respectively. 



5. Experimental results 

Figure [8] displays the dynamics of the amplitude determined from equation (J7J in the 
full range of A and B. The colour indicates sign and strength of the "force" f(B,A) 
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Figure 9. The multiquadric RBF approximation of dA/dt(A, B). The colour indicates 
this velocity in mm/s. The red dots correspond to the measured data. The blue crosses 
show the position of the centres of the radial basis functions. The black solid (dashed) 
line represents the stable (unstable) zero of the approximant. 



and confirms the bifurcation diagram suggested in figure HJ However, the data points 
are not only affected by noise, they are also irregularly distributed. The density is high 
at the stable solution branch, and on the contrary there are voids in the diagram, where 
no data is available at all. These voids belong to data sets which have been sorted out, 
because they converge to a localized solution (cf. § 14. 2p . 

In order to further interpret the data, a smooth approximation is helpful. The next 
section shows a way to construct such an approximation. Thereafter it is compared with 
analytical model equations (§ I5.2J) 

5.1. Approximation with radial basis functions 

We use a multiquadric radial ftasis /unction (RBF) network to compute a smooth 
sensible approximation to the amplitude equation (Powell 1990). A RBF network is a 
linear combination of shifted basis functions 



six 



E*= 



x 



p(x) 



(8) 



with the weights Xj, a basis function <ft, the arbitrarily chosen centres Cj and the low 
order multivariate polynomial p(x). We choose the multiquadric basis function 

<j>{x) = Vl + x 2 (9) 

and a linear polynomial p(x). The vector a? is a scaled combination of the coordinates 
of the phase space 



x = (B/5 b ,A/5a) 



(10) 
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5 B = 0.2mT 

8a — 1 mm. 

The scale variables Sa,b are chosen according to the empirical rule, that they should be 
approximately equal to the distance between centres in the corresponding direction. 

RBF networks with the basis function (Q are infinitely different iable and thus well- 
suited to the approximation of smooth functions. This kind of approximation has been 
proven useful in many applications. For an overview, see the review by Hardy (1990). 

Figure [9] displays the result from fitting the measured data from figure [8] with 
equation (jHJ). The centres Cj can be chosen almost arbitrarily; we have selected 50 
centres equally distributed in the measured region and additionally placed 68 centres 
near the stable solution branch. The location of these 118 centres is marked in figure [9] 
by the blue crosses. 

Only the weights Xj in equation (jSJ) are adjusted, while the position of the centres Cj 
are held constant. The problem of adjusting equation (JSJ) to our data is therefore reduced 
to a general linear least squares fitting problem with 118 free parameters. In order to 
avoid the problem of overfitting, we do not fit the data directly with equation (JHJ), but 
instead apply Tikhonov regularization (Neumaier 1998) 

J2\ S (x j )-f(x j )\ 2 + ^ X l = min! > ( n ) 

j k 

where f(xj) are the measured data and £ is the regularization parameter. This 
parameter controls the complexity of the model. A value near zero results in a standard 
least squares fit of the model to the data, while a larger value makes the resulting 
model more smooth. To find the optimal value for the regularization parameter, we 
apply standard 2- fold cross validation (Picard & Cook 1984), which yields an optimal 
regularization parameter of £ = 0.4. 

The zero from this fit gives an estimate for the stable and unstable solution of 
the amplitude equation and is indicated by the solid and dashed black line in figure ES 
respectively. For comparison, the original data are also shown with the red dots. The 
RBF reconstruction of the solution gives a very plausible result, which captures the 
essential features of the raw data, namely the imperfection, and provides an estimate 
for the unstable solution, which can only be approached to a certain limit by this method 
(see the voids in diagram EJ). This result of the RBF reconstruction has also been used 
to display the equilibrium in the figures 0] and [8] above. 

Such an approximation of the measured data by a RBF network has proven useful 
to get an unbiased estimate of the equilibrium and to fill in missing data. However, 
such a non-parametric regression technique cannot provide physically meaningful 
parameters (Green & Silverman 1994). In the next section, we compare the experimental 
data with nonlinear model equations for the bifurcation. 
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5.2. Comparison to nonlinear model equations 

A model equation for the imperfect, hysteretic bifurcation diagram like in figure H] is 
given by 

BA 

r — = £ A + ll (l + l2 e) A 2 - gA 3 + b, (12) 

where r denotes a time scale and e = (B 2 — B 2 )/ B 2 measures the distance from the 
bifurcation point (Reimann et al. 2005). For the diagram depicted in figure HI 71 must 
be positive, and the width of the hysteresis increases with the magnitude of 71 . In order 
to get stable patterns for e > 0, g must be positive, and the final amplitude decreases 
with the magnitude of g. The constant term b represents the imperfection of the system. 
For b > 0, the transcritical bifurcation at e = 0, A = dissolves into two saddle- node 
bifurcations, the distance of which is controlled by the magnitude of b. 

In the perfect case (6 = 0) and for 72 = equation ( fT2|) simplifies to equation (CQ). 
For 72 = 1 equation ( ITS]) comprises the amplitude equation (T5]), where the quadratic term 
is augmented by the factor (1 + e). From a formal point of view, one might ask, whether 
this additional factor could also incorporate another adjustable parameter, the blend 
parameter 72, resulting in the more general equation ffT2|) . Therefore, a comparison of 
the experimental data with this equation provides some hint, whether it is empirically 
advantageous to incorporate the e-dependence in the quadratic coefficient. 

The coefficients in the equations ([TJ) and ([2]) are in principle given by Bohlius 
et al. (2008) and Friedrichs & Engel (2001), respectively, as functions of the material 
parameters. However, these analyses are carried out for a linear magnetization curve 
M(H) = xH and, in the latter case, small values of x om y 5 an d therefore cannot be 
directly applied to our fluid. By treating r , B c , 71, 72, g, and b as adjustable parameters, 
these equations can nevertheless be fitted to the measured data. The coefficients are 
listed in table [2j 

Figure [TOl displays the result of fitting the model equations to the experimental data 
together with the neutral curve of the RBF approximation (white line) for reference. At 
a first glance, all models give reasonable approximations to the experimental data and 
the differences are minor. One such difference becomes apparent when looking at the 
equilibrium line (black solid line). The experimental data (red dots) evidently contain 
the final static amplitude in almost all cases, so the equilibrium should pass through, or 
very nearby the endpoints of the experimental data lines. This condition does not hold 
for the plain cubic equation in figure [TOT a). The width of the bistable regime ASh ys t 
suggested by the experimental data amounts to 0.29 mT, but the cubic equation provides 
only 0.18 mT (c.f. table [2]). The equation (T5]), which is augmented with (1 + s), gives 
a much better approximation of the equilibrium line (figure [T0l6). Here, the width of 
the hysteresis amounts to 0.25 mT, which is better but still smaller than the true value. 
Finally equation (fT2|) . which is a blend of the two former models, scores in between 
them in this test with Ai?h ys t = 0.24 mT (figure [TOb). The comparison of Ai?h ys t takes 
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Figure 10. Comparison of different models for the amplitude equation, (a) cubic 
equation (TTJ). (6) augmented equation @, and (c) blended equation (|T2)) with 72 = 
0.755. The white line is the neutral curve obtained from the RBF approximation, while 
the black line and the colour code comes from the corresponding model. 
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Table 2. The different models and their coefficients (columns 1-5). To elucidate their 
applicability to the experimental data we compare the width A_Bh yst of the bistable 
regime, the integral deviation H of the static solutions from the RBF approximation 
in the investigated range A £ [0.04 mm. . . 2.97 mm], and the residual sum of squares 
X 2 of all time dependent data (columns 6-8). For the amplitude equations, also the 
critical threshold B c is given in the last column. 



Model 


TO 

(s) 


7i 
(l/mm) 


72 


9 
(l/mm 2 ) 


b 
(mm) 


Ai?hyst 

(mT) 


10~ 6 


x 2 

(mm/s) 2 


B c 

(mT) 


Eq.fH 


0.4897 


0.0773 





0.0297 


0.0013 


0.1809 


7.58 


327.638 


11.27 


Eq. © 


0.5933 


0.1171 


1 


0.0407 


0.0024 


0.2468 


2.56 


326.834 


11.40 


Eq. mi) 


0.5770 


0.1111 


0.891 


0.0391 


0.0023 


0.2372 


2.85 


326.816 


11.39 


RBF 


N/A 


N/A 


N/A 


N/A 


N/A 


0.2834 


N/A 


278.12 


N/A 



only the two saddle-node points into account. To compare the whole equilibrium line, 
we provide an integral measure, namely the deviation in B on the static solution 

j4max 

K ■-- . '" ,.., / (B m v - B? cL4, (13) 



A ■ R 2 



J (B RBF -B) 2 dA, 



which is integrated over the whole amplitude range covered in the experiment. Here, 
-Brbf serves as a smooth approximation to the experimental data, and B c = 11.36 mT 
is the mean value of the critical inductions of the three model equations. 

For a list of K, see the 7th column of table |2j In this measure, the augmented 
equation (T5J) comes closer to the RBF by a factor of 3. Both inspected measures, Ai?h ys t 
and K demonstrate that the augmented equation (jSJ) provides a better approximation 
than the pure cubic equation (CQ). A description by the blended model ( fl~2l) yields a 
best fit blend parameter 7 2 = 0.891. This tends more to the augmented equation (j2j) 
than to the plain cubic equation ([T|). It is thus safe to conclude, that equation (j5j) 
provides a wider range of applicability compared to equation (JTJ, with the same number 
of adjustable parameters. 

Even though the equilibrium curve can be represented reasonably well by at least 
the augmented amplitude equation, the fitted models all deviate from the measured 
data off the equilibrium curve in two ways. At small inductions, the model equations 
underestimate the velocity of the pattern decay. Figure UTTa) depicts the amplitude 
equations for an undercritical induction of B — 10.80 mT. Here, the model equations all 
yield approximately the same result, but the difference to the measured data is roughly 
30%. At overcritical inductions (figure [TTTM. the model equations culminate around 
A = 2mm, while the RBF-approximation has its maximum around 1.5 mm, which 
coincides with the experimental observation. 

The quality of the theoretical descriptions becomes especially apparent when 
comparing the temporal evolution of the amplitude with numerical solutions of the 
amplitude equations, as shown in figure [12j These solutions are obtained by solving for 
A in equations (11121121) . with the first experimental data point as the initial condition. 
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Figure 11. Comparison of the amplitude equations for a constant induction B = 
10.80 mT (a) and B = 11.35 mT (6). The symbols display the original measured data, 
the dashed line is the RBF approximation, and the solid lines represent the parametric 
model equations ([TJ red), (0 blue) and (fT2l green). 



S 
S 




10 15 20 
Time t (s) 



30 



Figure 12. Comparison of the time evolution of the amplitude for the paths 2a and 
2b with the integrated model equations. The black dots display the original measured 
data, the lines represent the solution of equation ^ for different right hand sides, 
namely RBF (red dashed line) and the amplitude equations (TTJ red) , ([5J blue) and (fT2l 
green) . 



The solutions in the bistable regime (path 2a) describe the relaxation well, whereas the 
path 2b, observed at lower induction, differs notably from the integrated solution. This 
difference results from the underestimation of the decay rate far from the critical point, 
which was discussed above. Of course, the model equations are expansions about the 
critical point, and are not expected to be accurate far from the critical induction. In 
contrast, the temporal relaxation obtained from the RBF network fits the experimental 
data well for both paths. 
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A quantitative measure for the deviations of the model equations s(xj) from the 
experimental data f(Xj) is given by the sum of squared residuals 

x 2 = 5>(*.) -/(*;) I 2 - (14) 

3 

For the computation of x 2 , the dynamic evolution of the amplitudes in the whole range of 
the control parameter has been taken into account. The resulting values of x 2 are listed 
in the 8th column of table [2j The values obtained for the equations (fTf2][T2]) are larger 
than the one obtained for the RBF approximation. This indicates, that the dynamics 
away from the equilibrium curve is better described by the latter approximation. Note 
that in contrast to K, x 2 cannot discriminate between the different model equations 

(HE ED. 

For the overcritical induction, also the RBF approximation shows artefacts, as can 
be seen in figure [TTT 6). The measured data suggest a smooth shape of the amplitude 
equation, while the RBF approximation shows a bump near A = 2... 2. 5 mm. This 
irregularity stems from the double equilibrium point, which can be seen in figure [6] 
(c.f. traces 2a and 3a) and figure [71 These double equilibrium points represent similar, 
but different stable patterns and exist also in the neighbourhood of the data shown in 
figure [TTT6). They influence the RBF network, because it is a global approximation 
scheme. 

6. Discussion and Outlook 

Using a highly viscous ferrofluid, the pattern evolution in the succession of the 
Rosensweig instability can be slowed down to the order of minutes. Therefore, it 
is possible to conveniently observe the dynamics using a two-dimensional imaging 
technique based on the absorption of X-rays. We have explored the full phase space 
in the vicinity of a hysteresis loop formed by a transcritical bifurcation and a saddle- 
node. From the measured data we have reconstructed a nonlinear equation of motion. A 
smooth approximation of the noisy experimental data with radial basis functions (RBF) 
provides an estimate for missing data and for the static amplitude. Amplitude equations 
for the instability were successfully fitted to the experimental data. The estimate of the 
equilibrium curve based on the cubic amplitude equation, as derived by Bohlius et al. 
(2008), can be significantly improved by augmenting the quadratic coefficient with the 
factor (1 + e). This factor can be extracted from the free energy functional of Friedrichs 
& Engel (2001) and reduces the deviation K between the equilibrium curve and the RBF 
reference by a factor of three. 

However, a comparison of the absolute magnitude of both amplitude equations with 
the experimental data at an induction 5 % below the critical induction B c reveals, that 
neither amplitude equation can exactly describe all of the recorded dynamics. At least 
the agreement between the experiment and the model equations is within 30 %. This 
is surprising, when considering that these equations are low order expansions about the 
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critical threshold B c , while the investigated range of about 7.5 % of B c even includes 
another critical point at a saddle node. 

So far, the experimentally determined expansion coefficients can not sensibly be 
compared with the ones derived by theory, because the latter does not include a nonlinear 
magnetisation curve, which has a major influence on the instability. For example, for 
our fluid the critical threshold B c is shifted by about 10%. The change of the expansion 
coefficients is unknown. Its calculation remains a challenge for theory. 

While we have extracted a nonlinear equation of motion for the specific example of 
the Rosensweig instability, it seems worth mentioning that the method is applicable to 
a wide set of overdamped pattern forming systems. A necessary precondition is that the 
control parameter can be switched much faster than the characteristic relaxation time. 

In this paper only space filling homogeneous patterns have been considered for 
the method. However, the real system has additional degrees of freedom, which are 
manifested in the formation of localized states as observed previously by Richter & 
Barashenkov (2005). Whereas in the latter work a local disturbance of the magnetic 
induction was needed to ignite localized states, here localized spikes and hexagonal 
patches emerge spontaneously when the measurement protocol pushes the system into 
the neighbourhood of the unstable branch. While a nuisance for the particular task 
of this paper, this multistability of hexagonal patches is an interesting phenomenon in 
itself, and may be seen in connection with homoclinic snaking (see e.g. Lloyd et al. 2008). 
However, for a clear picture further measurements are needed. 
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